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Abstract 

The evolution equation for q q production introduced by Marchesini and Mueller 
posed some intriguing mathematical puzzles, both numerical and analytic. I give 
a detailed account of the numerical approach which led eventually to the exact 
solution. While part of the work was in fact along a wrong track, it turns out that 
some of the techniques involved are interesting in their own and applicable to many 
other problems, i.e. to the numerical study of Ricci flows. 
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1. Introduction 

Marchesini and Mueller ^ |2| introduced the following equation for the 
evolution of q q evolution in QCD 

du(r,g) _ f l u(T,r)£)/ri -u(t,Q ^ | f l u(t,£,/i]) - u(t,£) ^ ^ ^ 



dr J 1 - 77 Jz 1-rj 

where the unknown function u(r, £) must vanish at £ = to ensure the 
convergence of the integrals involved. If we put = £ip(£) (boundary 
conditions are then taken into account if ip is bounded or at least does not 
grow too rapidly at £ — ► 0) 

= (K4,)(t,Z)=£ -^ ^( T ,T,)-mm(l,ar } MT,0) (1-2) 

I shall refer to K as the MM operator. We can easily discretize K on a lattice 
£n = no, a = 2/{N + 1) 

w) t = y ti^MMi . (i.3) 

^-^i^ \i ';— j\ 

Any trace of a disappears from the discrete equation which is a sign of the 
scale invariance of the original equation: under £ — > A£ only the endpoint 
changes, hence the result is insensitive to its actual value. The spectrum 
can be estimated numerically. By Richardson extrapolation from dimension 
32, 64, ... , 1024 one gets for the spectrum of K 

E = 

2.4990 

1.7993 

0.9604 

0.2179 
-0.3663 



K has a negative spectrum except for a few positive eigenvalues, the 
largest one dominates the evolution (all others are damped away). 

Notice that if (accidentally, by mistake!), you ignore the "min(l, £/r/)" 
factor in Eq. (|1.2jl the spectrum comes out very simple and, surprisingly 
enough, independent from N to all available figures: 
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What is the secret behind these numbers? Taking the differences we get 

2.0000 
1.0000 
0.6667 
0.5000 
0.4000 
0.3333 
0.2857 
0.2500 

an easily recognizable sequence. Indeed the eigenvalues are given precisely 
by twice the "harmonic numbers" {t) n = Y^j=i J 1 I " > 0} and f)o = 0. 
This fact is actually an exact property of the modified integral equation, 
both in its discretized form and on the continuum (a result which goes 
back to the sixties (SI, see Appendix), the eigenvectors being Tchebyshev 
discrete polynomials |IJ Ej which converge to Legendre polynomials in the 
limit N — > co. 

2. Perturbation theory 

Following the hint of the previous Section, let us represent K as the sum 
of two terms and treat the problem by perturbation theory. 

W)(0 = {KoMQ - log(£) m (2-1) 

where 

{Kom)= f\,<!^z*^ (2 . 2) 

Jo 

Kq is exactly diagonalizable, with eigenfunctions the Legendre polynomials 
P n (2£ — 1) and eigenvalues proportional to the harmonic numbers h n = 
Y^=i i _1 ( see Appendix A). Second order perturbation theory gives for the 
ground state Eq ~ 1.44754 hence convergence appears to be rather slow. 
The usual methods to get high order coefficients are not applicable here, 
since the matrix (P n \ log(£)|P m ) is full. 

One can do better with a purely numerical approach as we discuss in the 
next section (the coefficients, we shall see, decrease only as l/nlog(ra), which 
would require high orders in p.t. to get a meaningful result). 

3. Evolution 

I recyvled an old program which was used to study the renormalization 
group equation of the non-linear sigma model (the "sausage" [El III IE])- The 
equation is now rather popular in the mathematical literature as the Ricci 
flow. The idea is to split the evolution of 



(3.1) 
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into two steps 



^(A + r,0 = Vw(6 ~ r log(e) V(A, 0) 



(3.2) 



The first step is accomplished by going to the representation in terms of 
Legendre polynomials (tp = ^ ip n -P n (2£ — 1)) where Kq is diagonal. Coming 
back to the ^-representation one executes the second step. The program is 
implemented in matlab. 




Figure 1 . Evolution from £ = 1 



En passant one can study the spectrum of K within the same program. 
We get the result of Tab OH where the last line is obtained by extrapolating 
in the variable ra" 1 / 4 , which appears at first sight as an approximate scaling 
law (but see later on). This value should be compared to the approximate 
saddle point value 4 In 2 ~ 2.77. The strong dependence on the grid size is 
not surprising, since we have to deal with a singular scale invariant integral 
operator. The similar operator considered by Tuck is not scale invariant 
and its cutoff dependence is much flatter (i.e. better!). 
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Table 1. K ground state by the direct Legendre transform. 



4. Alpert-Rokhlin's Fast Legendre Transform 

Due to the very slow convergence toward N— > oo it is desirable to be able 
to calculate the spectrum with a high number of collocation points. This is 
totally unfeasible with the direct method. The calculation with n = 4096 
required a work space of 1/2 GByte and going further was not possible on 
available workstations. The way out is to apply some sparse matrix com- 
putational tool which should be able to save memory and time. It was 
shown by Alpert and Rokhlin .1 j)] that it is possible to transform from a Le- 
gendre expansion Yln<N °n Pn{x) to a Tchebyshev expansion J2n<N &n T n {x) 
in 0{N) time, even if the amount of memory required may be rather large 
(at least 200 N words). Since Tchebyshev polynomials of the first kind 
are just trigonometric functions in disguise, the Legendre transform is re- 
duced to a combination of Alpert-Rokhlin's transform (art) and cosine- 
Fourier-transform. Using Alpert's implementation of art 2 combined with 
fftw in mode REDFT10/01 (see Ref.|llj) 3 we realized a code essentially 
equivalent to the previous one but allowing for high dimensional matrix 
representation of the operator. The ground state has been computed for 
N = 2 k ,k = 6,7,... ,18 giving the result of Tab. 2. The difference from 
the previous calculation is due to a different choice of discretization grid 
(Gaussian integration points, i.e. the roots of -Pjv(x), in the former case, 
Tchebyshev points, uniformly spaced in acos(x), in the latter). Notice that 
the results of the fast method anticipate those of the direct method, that is 
the "fast" result at N is close to the "direct" result at 2N. In a sense the 
formal dimension of the real fftw (2N) is the "true" dimension. 

For the technically-oriented reader we report the approximate timings of 
the two algorithms in Appendix B. 



B. Alpert very kindly provided us with his Fortran code, 
'http:/ /www. fftw. org/#documentation 
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The extrapolation at A r — > oo , assuming a power law scaling as before, 
seems consistent, giving 2.6733 and 2.6692 (linear and quadratic fit respec- 
tively) with the first method, 2.6661 and 2.6704 with the second. We would 
conclude that the saddle point estimate 41og(2) is correct within 4%. 
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Table 2. The ground state from Direct Legendre Transform 
(DLT) and from the fast algorithm (art+fftw) 



A totally different result is however hiding behind these figures. It must 
be realized that the crucial point is to identify the correct N dependence, 
since this is going to make a big difference in the extrapolation at N — > 
oo. A careful analysis shows that a logarithmic scaling law is much 
more accurate than a power law. Looking for a fit of the kind E(N) ~ 
E(oo) — Ci/log(A r /A) — Cil log(iV/A) 2 we get a very good interpolation 
(the deviation is uniformly less than 1 part in 10 4 ) and the value at N = oo 
is compatible with 41og(2) (within the same accuracy). According to this 
idea we should conclude that, surprisingly enough, the saddle point value is 
actually exact (see Figg. 3,4). In the case studied by Tuck we find a much 
steeper scaling law of the kind 

E(N) = Ei(oo) + CN~ 2 log(N)- 1 

as shown in Fig. 5. 

It has been realized [2] that the picture is simply due to the different 
character of the spectrum: continuous for MM and discrete for Tuck's op- 
erators. This fact is made absolutely transparent by adopting a different 
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Figure 2. Continuum limit for the ground state with a 
power scaling law. 



representation of the operator K: 

{K<j))(x) = {K Q (j>){x) +21og(l + e^ /2 )0(x) -41og2 0(s) 

{K <f>)(x)= . — — -dy. 

Jo 2smh| 5 (x - y)\ 

The integral operator Kq is almost local and it is not very different from a 
kinetic term. If we consider a wave-function with support in a region far 
from the origin, the operator reduces to 

J — OO ^ olLLLL | 2 V X a ) | 

which is diagonal in Fourier space with eigenvalue — the subtracted 
Lipatov function, explicitly given by x(p) = 2(-0(l) — K(V'(2 1 + w))) ~ 4 log 2, 
ip being the logarithmic derivative of the V function. Kq is well-known, not 
necessarily in this form, as the BFKL operator |12j . 

It has been realized that the representation introduced here is also more 
convenient to allow a numerical study of the evolution in the case of MM, 
while this is not the case for Tuck's equation. Essentially the dominant 
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FlGURE 3. Continuum limit with the logarithmic scaling law. 
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Figure 4. Extrapolation in log-scale shows the emergence 
of a continuous spectrum 



eigenvalue 4 log 2 is already built-in, while in the representation of SecEJthis 
value can only be obtained by extrapolating at very large matrix dimensions 
(see Tab. ij). 
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Figure 5. Continuum limit with the modified scaling law, 
Tuck's case. 



To make the difference between mm and bfkl operators more explicit, 
it would be desirable to be able to apply the method of images (which is 
usually employed with local differential operators) to get rid of the boundary. 
However no simple boundary condition seems to be appropriate. Actually by 
solving the eigenvalue equation by standard linear algorithms (matlab's eig 
routine) one finds that the eigenvectors are essentially shifted trigonometric 
functions, i.e., far from the boundary, 4>k(x) ~ sin(/cx + 5(h)) . By switching 
V{x) on and off, we can easily check that the behavior at x = is strongly 
influenced by V(x). 

The phase shift 5(k) is particularly interesting. For example, the asymp- 
totic behaviour of </>(A, x) at large A is strongly influenced by it. This fact 
is well-known in the theory of potential scattering in quantum mechanics. 
While the general setup here is quite different, nonetheless there are remark- 
able analogies which give useful guidelines. For example the vanishing of 5 
at k = is a signal of the absence of bound states (Levinson's theorem), 
were we able to extend the theorem to this context. Details can be found in 

0- 

The "unbounded" representation helps in understanding what really goes 
wrong with the initial approach based on the Legendre basis. Introducing 
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a finite box of side L (0 ^ x ^ L) the energy spectrum is discretized and 
at low energy it is given by oc (nir/L) 2 . In the Legendre expansion of the 
previous section all Gaussian points are confined to x < L = 2 log(iV). This 
fact explains the logarithmic scaling law depicted in Fig. 3. Also, since a 
good description of the evolution at large r requires L 3> 100, this cannot 
be explored through the Legendre expansion. 

5. Further developments 

Recent developments pushed our understanding of the problem to a higher 
level. A precise characterization of the time-dependent solution of MM equa- 
tion was developed by a perturbative technique which can be pushed to all 
orders and allows for a full resummation An exact form for the phase 
shift and the continuum eigenfunctions has been derived. A rigorous proof 
on purely algebraic grounds, thus avoiding a delicate problem of resumma- 
tion, has been later found, thanks to an idea of V. A. Fateev |13| . 

There exists another representation of the integral operator which avoids 
the presence of a boundary. Thanks to the intrinsic scale invariance, the 
equation can be remapped on the whole of R by setting exp(x) = exp(x') + 1, 
which leaves the kernel invariant and only modifies the potential. In this 
representation we may apply a spectral algorithm to the evolution equa- 
tion simply based on Fourier transform, more economic than the combined 
art+fftw. This will be left homework. 

6. Conclusions and outlook 

The integral equation introduced by Marchesini and Mueller is deeply 
related to another problem in mathematical physics studied by E. Tuck 
fourty years ago. The connection to Tuck's equation was used to analyze 
MM operator's spectral properties by an efficient (sparse) matrix computa- 
tion, based on Alpert-Rokhlin transform, fftw and arpack. This analysis 
suggests that the spectrum is continuous with endpoint 4 log 2. A second 
representation of the integral operator makes the spectral properties more 
transparent and lends itself to an easier algorithmic implementation which 
allows to evaluate the evolution at large A. 

The application of art to the renormalization group equation for 0(3) a- 
model may be useful to achieve greater accuracy than allowed from the direct 
transform 8 . More generally, the application of a full group theoretical 
transform without axial symmetry, will make it possible to explore the 0(3) 
Ricci flow in full detail. 
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Appendix A 

Here is a proof for the discrete form of Tuck's operator. For the original 
problem see El • 
Let 

I tA n ) \ V^' V i ~ Vi fm\ 
(K > v) i = ^ -j— — . (6.1) 

Let's apply the matrix Kq N ^ to the vector vf = i e . We have 

= (-y- 1 +y N 

^ ^j=o ^j=i+i J ^k=o (6.2) 

^k=o \ k + i v 7 

where we used the Euler-MacLaurin summation formula for ^ . j k , and the 

constants {c n } are calculable but unnecessary This proves that KW v^' is 
contained in the linear span of [v^°\ v^]. Since Kq(N) is symmetric, 

it is diagonable, its eigenvectors are orthogonal, hence they are given by the 
orthogonal discrete polynomials with respect to the uniform weight on the 
set [0, 1, 2, . . . ,N]. The eigenvalues can be read off the coefficient of in 

the expansion of Kq N ^ . The explicit form of the eigenvectors is given by 
Tchebyshev polynomials of a discrete variable 1,4.. 



Appendix B 

We give here some technical details about the algorithms which we have 
applied in the paper. The direct method consists in building the table of 
Legendre polynomials {P n (xj ), n = 0, 1, . . . , N-l} at the Gauss points, i.e. 
at the roots of Pn(x). The technique, exploiting the recurrence relation of 
orthogonal polynomials, is due to Golub and Welsch (see Ref.|15j). To find 
the spectrum of K we simply define Kq to be diagonal in the basis {P n (x)} 
with eigenvalues —1\\ n . The matlab routine "eig" is then invoked. The sin- 
gularity of the logarithm at the boundary is avoided because the zeros of 
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the polynomials are all internal at the interval [—1,1]. The matrix repre- 
sentation of the free part Kq is exact, since the Gauss quadrature formula 
is exact on polynomials of low degree. In finite precision arithmetic Kq is 
affected by the accumulation of truncation errors, yielding an error of order 
1CP 13 on its spectrum, which is rather irrelevant. The method is presently 
feasible for dimension less than 4000 and it has the advantage that N can 
be any integer, not necessarily a power of 2. 

The method based on ART makes use of the expansion on Tchebyshev's 
polynomials of the first kind T n (x) = cos(n arccos(x)). Again the Gauss 
points are interior at the interval and the singularity is avoided. The real 
DFT of kind REDFT10 precisely makes use of this grid of points. Even 
if the Gauss-Tchebyshev integration is exact for polynomials of low degree, 
still a problem arises, namely that Kq is a symmetric operator with respect 
to the Lebesgue measure whereas Tchebyshev's polynomials are orthogonal 
with respect to a different measure. It turns out that to restore symmetry 
we have to deal with Kq = (1 — a:) 1 / 4 Kq(1 — x) -1 / 4 , hence we are outside 
any polynomial subspace and this introduces a systematic error making Kq 
only approximately symmetric. This has been checked after realizing that 
the spectrum of Kq considered as a symmetric operator contains substantial 
error, up to 10%. The strategy we adopt is therefore to relax the symmetry 
condition on Kq and compute the spectrum with a version of the Arnoldi 
algorithm which allows to deal with non-symmetric operators provided by 
arpack[TH|. The combination ART+FFTW+ARPACK(this latter offers the 
routine znaupd which applies to general non-symmetric complex matrices) 
turns out to be again in the game with an accuracy on the spectrum of Kq 
comparable if not superior to the direct method. For example at N = 64 
we find 

E=[ 

1.00047386511033e-18 
2 . 00000000000000e+00 
2 . 99999999999999e+00 
3 . 66666666666666e+00 
4. 16666666666666e+00 
4 . 56666666666665e+00 
4 . 90000000000001e+00 
5.18571428571429e+00 

]; 

and the reader can verify by herself that the error is only at the last decimal 
place. A further (marginal) improvement will be achieved using the real- 
non-symmetric routine dnaupd; we use the complex version because it was 
already implemented as a C++ header by Scott Shaw. 

Let us now comment upon performance. The first methods grows in time 
and memory rather quickly (0(N 2 ) in memory and 0(iV 3 ) in execution time. 
The fast one is much less memory greedy. Notice that in the largest case 
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examined, N = 2 , the program requires slightly less than half a Gigabyte 
of memory, half of which simply to allocate art's working array (200./V 
words). By contrast the direct method requires ~ 450 MBy at N = 4096 
and it would grow to the order of 6 TBytes at N = 2 18 while the execution 
would require 50 years at the present cpu speed. 

Tables 3,4 report the execution times on a pentium III with clock at 1.13 
GHz and on a Xeon with clock at 2.8 GHz, respectively, using Matlab v. 6. 5. 
Execution times for the "fast" algorithm are inclusive of the preconditioning 
(we select the initial vector by executing a number of Trotter steps). As 
it is rather clear from the table, the execution time grows as expected as 
0(Nlog(N), with some expected deviations when the system switches to 
virtual memory. 
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Figure 6. Execution times for the Direct Legendre trans- 
form and Alpert-Rokhlin transform. 
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